UNIVERSITY  OF  MARYLAND 
COMPUTER  SCIENCE  CENTER 

COLLEGE  PARK,  MARYLAND 
20742 


for  jMfcU*  nlNH  t 


* 


} 


j 

1 


i 


am  ftmcs  office  of  scientific  iamr * 

NOTICE  OF  TRANSMITTAL  TO  MO  ^ 

This  technical  report  has  been  reviewed  4. 
approved  for  public  release  LAM  UR  190-lM  17b) 
Distribution  is  unllalted.  '  D' 

A.  0.  BLOSfi 

Teohnloal  Information  Officer 


7^ 


TR-909  /  - 
TaF  GSR-')  7  -"327 1 1 


/  STOCHASTIC  MODELS - 

i  - ■'  FOR  pLOSED  BOUNDARY  ANALYSIS: 

/  PART  I  r  REPRESENTATION  AND  RECONSTRUCTION 

■f~r.Y>  R.  L./Kashyap 
'\'J.  ft.ycheiiappa^'  + 


jU  s 


— 


/  /  v/V/)/w 


ABSTRACT 


This  paper  deals  with  the  analysis  of  closed  boundaries  of 
arbitrary  shape  in  a  plane.  Specifically,  it  is  concerned  with 
the  problems  of  representation  and  reconstruction.  We  first  set 
up  a  one  to  one  correspondence  between  the  given  closed  boundary 
and  a  univariate  or  multivariate  sequence  of  real  numbers.  Uni¬ 
variate  or  multivariate  circular  autoregressive  models  are  sug¬ 
gested  for  the  representation  of  the  sequence  of  numbers  derived 
from  the  closed  boundary.  The  stochastic  model  representing  the 
closed  boundary  is  invariant  to  transformations  of  the  boundary 
such  as  scaling,  rotation  and  choice  of  the  starting  point. 

Methods\  for  estimating  the  unknown  parameters  of  the  model  are 
giventand  a  decision  rule  for  choosing  the  appropriate  order  of 
the  model  is  included.  Constraints  on  the  estimates  are  derived 
so  that  the  estimates  are  invariant  to  transformations  of  the 
boundaries.  The  specific  stochastic  model  used  enables  us  to 
reconstruct  a  closed  boundary  with  less  computational  effort  using 
FFT  algorithms.  Results  of  simulations  are  included  and  applica¬ 
tions  to  contour  coding  are  discussed.  In  a  subsequent  paper  we 
will  consider  the  classification  problem.  <- 
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1 .  Introduction 


The  problems  of  closed  boundary  analysis  and  discrimination 
are  familiar  in  pattern  recognition.  They  arise  in  a  variety  of 
situations  such  as  contour  coding  [1,2],  classification  of  chromo¬ 
somes  [3] ,  interpretation  of  chest  x-rays  [4] ,  scene  analysis 
[5],  and  the  detection  of  parts  in  a  bin  in  robotics  [6].  In 
this  first  part  of  our  study  of  the  applicability  of  stochastic 
models  to  such  problems,  we  consider  the  analysis  of  boundaries. 

In  a  subsequent  paper,  classification  of  a  closed  boundary  as 
one  of  several  mutually  exclusive  boundaries  will  be  considered. 
Specifically,  in  this  paper  we  consider  the  problems  associated 
with  the  representation  and  reconstruction  of  closed  boundaries. 

The  statistical  approach  to  shape  analysis  has  so  far  been 
concerned  with  the  Fourier  analysis  of  some  function  derived  from 
the  boundary  [7-10].  This  is  essentially  a  deterministic  approach 
Moreover,  in  practical  applications,  the  Fourier  series  expansion 
is  truncated  to  a  certain  number  of  terms,  leading  to  an  inexact 
representation  of  the  data.  Further,  not  all  Fourier  descriptions 
yield  coefficients  which  are  insensitive  to  rotation,  scaling 
and  variability  in  the  starting  point  in  the  tracing  of  the 
boundary  [7]  and  the  coefficients  do  not  always  regenerate  a 
closed  boundary  [8]. 

In  this  paper  we  suggest  a  stochastic  model  of  a  special 
type,  the  so-called  Circular  AutoRegressive  (CAR)  models,  to 


describe  the  data  from  the  boundary.  The  given  boundary  is 
approximated  by  a  series  of  straight  line  segments  and  the  two- 
variable  data  formed  by  the  (x,y)  co-ordinates  of  the  end  points 
of  the  line  segments  denoted  as  { Xj^  (i)  ,x2  (i) } ,  i  =  1,2,...,N, 
are  used  for  the  representation.  One  of  the  characteristics  of 
the  data  is  its  circularity,  i.e.,  x^ (t+N)  =  x^ft),  x2(t+N)  = 
x2(t).  For  the  class  of  objects  in  which  each  of  N  radii  vectors 
from  the  centroid  of  the  boundary  intersect  the  boundary  at  only 
one  point,  a  simple  version  of  the  two-dimensional  representation 
is  used,  i.e.,  a  one-dimensional  series  of  real  numbers  is 
formed  by  measuring  the  lengths  of  successive  radii  vectors  that 
are  angularly  equispaced.  The  observations  {r (1) , . . . ,r (N) }  again 
possess  the  circularity  property. 

The  observation  set,  be  it  one -dimensional  or  two-dimensional, 
is  represented  by  CAR  models  characterized  by  a  set  of  unknown 
parameters.  By  making  suitable  assumptions  about  the  noise  driving 
the  model,  expressions  are  obtained  for  the  maximum  likelihood 
(m.l.)  estimates  of  the  parameters  of  the  model.  In  actual 
applications,  the  so-called  least  square  estimates  are  used  which 
are  an  approximation  to  the  m.l.  estimates.  The  conditions  for 
these  estimates  to  be  invariant  to  scaling,  translation  and  start¬ 
ing  point  on  the  boundary  are  derived.  It  also  turns  out  that 
these  estimates  are  the  sufficient  statistics  for  the  unknown 
parameters.  We  also  consider  the  reconstruction  problem  and 
present  algorithms  for  generating  a  closed  contour.  The  CAR 
models  yield  the  reconstruction  of  a  closed  curve  in  O(NlogN) 


operations.  The  data  compression  aspects  of  the  model  are  con¬ 
sidered  with  specific  applications  such  as  the  archival  storage 
of  weather  maps  [22]  and  contours. 

The  organization  of  the  paper  is  as  follows:  Section  2 
describes  the  different  ways  of  obtaining  suitable  representations 
from  a  closed  boundary.  The  general  representation  suggested  is 
two-dimensional  and  is  applicable  to  all  closed  boundaries  in  a 
plane.  A  one-dimensional  representation  is  used  for  the  class  of 
objects  in  which  each  of  the  N  radii  vectors  intersect  the  boun¬ 
dary  at  only  one  point.  In  later  sections  all  the  analytical 
results  are  obtained  for  the  case  of  a  one -dimensional  model, 
which  can  be  easily  generalized  to  two-dimensional  representations. 
The  mathematical  model  is  introduced  in  Section  3  and  the  proper¬ 
ties  possessed  by  the  estimates  of  the  parameters  characterizing 
the  model  are  described  in  Section  4 .  A  decision  rule  to  choose 
appropriate  neighbors  is  also  given.  In  Section  5  the  coding  of 
closed  boundaries  is  considered.  Results  of  the  simulation  are 


given  in  Section  6.  Discussion  is  given  in  Section  7.  Main  re¬ 


sults  are  stated  as  theorems , 
appendix . 


the  proofs  being  given  in  the 


Representation  of  Closed  Boundaries 

We  will  represent  the  given  closed  boundary  by  a  finite 
sequence  of  real  numbers,  the  so-called  time  series.  We  assume 
that  the  boundary  has  no  crossovers.  We  first  approximate  the 
boundary  by  a  polygon  of  N  sides.  The  time  series  exactly  repre¬ 
sents  the  polygon,  i.e.,  given  the  time  series,  we  can  reconstruct 
the  polygon,  and  vice  versa.  We  can  choose  the  integer  N  large 
enough  to  obtain  the  desired  accuracy  of  approximation.  We  em¬ 
phasize  that  the  time  series  is  constructed  such  that  the  corre¬ 
sponding  shape  analysis  and  classification  rules  are  invariant  to 
transformations  of  the  boundary  such  as  translation,  scaling, 
rotation  and  the  starting  point  in  the  tracing  of  the  boundary. 

We  now  describe  several  methods  of  obtaining  the  time  series  from 
the  boundary. 

Method  1  (Two-Dimensional) :  Consider  a  Cartesian  co-ordinate 
system  such  that  the  origin  is  at  the  centroid  of  the  object. 

Let  (xjd),  x2(i)},  i  =  1,2,...,N,  represent  the  co-ordinates  of 
a  large  number  (say)  &  of  points  on  the  boundary  as  in  Fig.  1. 

The  required  representation  is  { (x^ (i) ,x2 (i) )» i  =  1,2,...,N}. 

Method  2  (One-Dimensional)  :  Let { (x^ (i) ,x2 (i) ) ,  i  =  1,2,  ...,K} 
be  the  coordinates  of  a  large  number,  say  N,  of  points  on  the 
boundary.  Let 

N  N 

x.(0)  -  (1/N)  E  x. ( i) ,  x,(0)  =  (1/N)  I  x,(i) 
i*l  1  *  i=l  * 


and  let  (x^ (0) ,x2 (0) )  be  denoted  as  the  origin  0.  Let  A^  be 
any  point  on  the  boundary.  Mark  off  (N  -  1)  points  A2,A3,....A 
such  that  the  angles  between  0A^  and  0A^+^  are  all  the  same, 
viz.,  (2ir/N)  radians.  This  is  shown  in  Fig.  2.  Let  the 
distance  0A^  =  r(i),  i  =  1,2,...,N.  The  required  polygon  is 
A.  ,  A-  , .  . .  ,A  .  and  the  time  series  is  {r  (1)  ,r  (2)  ,  . . .  ,r  (N)  }  . 

1  Z  N 

Method  3  (One-Dimensional) :  We  can  construct  the  points  A^ , . . . , Aj 
as  in  Method  2.  Let  the  distance  between  A^  and  A^+^  be  b(i),  i  = 
1,2,..., N-l.  b(N)  is  the  distance  between  AN  and  A^ .  The  re¬ 
quired  time  series  is  {b (1) ,b (2) , . . . ,b (N) } . 

Method  4  (One-Dimensional) :  We  fit  a  polygon  of  N  sides  to  the 
boundary  such  that  all  the  sides  have  the  same  length.  Call  the 
polygon  AiA2*,,ANAi'  Let  <|k  be  the  angle  between  A^A^+1  and  a 
reference  line.  Then  {(Jk,  i  =  1,2,...,N}  is  the  required  time 
series . 

Note  that  methods  2 ,  3  and  4  can  be  used  to  represent  what  we 
shall  call  Type  1  objects.  The  class  of  Type  1  objects  include 
non-convex  objects  as  in  Figure  2. 

Definition  (Type  1) : 

In  any  type  1  object  the  N  radii  vectors  from  the  centroid 


of  the  object  intersect  the  boundary  at  only  one  point. 


Another  possible  representation  for  closed  boundaries  is 
the  tangent  angle  versus  arc  length  [8] .  One  disadvantage  of 


this  representation  is  that  the  angle  function  is  more  sensitive 
to  the  noise  inherent  in  a  fuzzy  boundary.  The  angle  function  is 
related  to  the  derivative  Of  the  coordinate  function,  and  hence 
small  variations  in  the  co-ordinate  values  of  the  boundary  points 
can  result  in  large  variations  in  the  direction  of  the  tangent 
vector. 


3.  Mathematical  Modeling  of  Boundaries 

Henceforth,  we  consider  the  analysis  with  respect  to  the 
one-dimensional  data  {r (1) , . . . ,r (N) }  derived  by  using  Method  2. 
These  results  are  generalizable  to  the  two-dimensional  represen¬ 
tations  of  the  non-Type  1  objects. 

We  are  given  a  one-dimensional  sequence  (r (1) ,r  (2) , . . . ,r (N) } 
such  that 

r(N+k)  =  r(k),  k  =  1,2,...,N 

r(-k)  =  r (N-k+1) ,  k  =  1,2,. . . ,N  (3.1) 

We  will  fit  a  particular  type  of  stochastic  process  called  a 
Circular  AutoRegressive  (CAR)  Process  as  in  (3.2): 


m 


r  (t)  =  oi+  E  6  r(t— t.)  =  /8u)(t),  t  =  1,...,N  (3.2) 

3=1  3  3 


where  t  ^ ,  j  =  1,2, . . . ,m  are  positive  integers. 


u)  (N+t )  =  <D(t),  t  =  1,2,. ..,N, 


E  (u)  (i) )  =  0,  E  (u>  (i )  u)  ( j  ) )  =  S'5  —  *  l^ijjs-N 


6.  .  /  0  if  i  =  j 
=  0  otherwise, 


<*>(•)  is  independent  and  normal  (0,8)  and  {o,6^,...,0  ,8} 


are  unknown. 

A  similar  stochastic  model  is  used  for  both  co-ordinates  of 


the  two-dimensional  representation,  i.e., 

m 


x^t)  = 


a.+  E  0.  .x.  (t-t.  .  )+/fTu>.  (t)  ,  i  =  1,2 
j=i  13  1  13  11 


4 


a 


4 


where  x^(t+N)  =  x^(t),  i  =  1,2  and  and^t*))  are  two  indepen¬ 

dent  normal  i.i.d.  noise  sequences.  This  simple  extension  of  (1.2] 
will  be  used  throughout  for  non-Type  1  objects.  By  introducing 
interdependence  between  x^(*)  and  X2 ( * )  more  general  models  can 
be  considered. 

Because  of  the  circularity  assumption,  given  u  (1)  , .  .  .  ,u>  (N)  , 
we  can  obtain  r (1) , . . . ,r (N)  by  solving  the  system  of  N  equations 
in  (m+1)  variables  given  in  (3.3): 


•°1  1 
■°2  -01 


-02  -0! 
9  3  -°2 
■°4  -03 


(3.3) 


0  . . .  -e„  -e  . 

m  m-1 


or  equivalently, 

B  (0  )  R  =  A  +  /3W 

where  R  =  [r  (1) , . . . ,r (N) ]T 

T 

A  =  (a,  ...»  a]  ,  Nxl  vector 
W  =  [u)(l)  , .  .  .  ,w  (N)  ]  T 

B ( 0 ) =  circulant  [1,0,0, ... ,-0,  -0m_^,-0^]. 

B (0 )  is  a  circulant  matrix  and  it  can  be  generated  from  its 
first  row  given  above.  The  probability  density  of  the  set  of 
observations  is  given  in  (3.4): 


p  (r  (1)  , . . .  ,r  (N)  |  a,  8 ,6) 

N 

=  det  |  B  ( 0 )  |  [1/2tt0]N/2  exp  {(-1/20)  Z  (r(t)  -  a  -  0T2T(  t-1 )  )  2 }  (3.4) 

t=l 

where  0  =  (0  , . .  .  ,  0  ) T 

~  1  m 

z(t-l)  =  (r (t-t^) ,  r (t-t2) , . . . ,r (t-t^) )T 
and  | det  B  (0)|  is  the  Jacobian  of  the  transformation  from 
r (1 ),..., r (N)  to  w  (1)  ,  .  .  .  ,u)  (N)  and  has  the  following  expression 
[11]  in  view  of  the  circulant  structure  of  B(0): 


I  det  B  <  6 )  1  A  |  B  (0)  |  =  IT  (1  -  )  (3.5) 

k=0  ~ 

where  the  m-vector  is  as  follows: 

=  column  [exp  (-/^127Tkti/N)  ,  i  =  l,2,...,m] 

In  our  formulation,  the  shape  of  the  boundary  is  controlled 
by  the  coefficients  (a,  0,0).  The  random  sequence  w(l),to(2),..., 
oi  (N)  represents  the  effect  of  noise  on  the  boundary.  For  instance, 
given  N  i.i.d.  random  variables  u>  (1)  ,  . .  .  ,u  (N)  from  a  Gaussian  or 
Gamma  distribution,  we  can  account  for  a  class  of  shapes  which 
are  minor  variations  of  the  basic  shape,  the  so-called  local 
variations  of  the  shape.  On  the  other  hand,  global  transforma¬ 
tions  like  scaling  and  rotation  of  the  given  shape  are  controlled 
by  the  parameters  (a, 0,6).  Note  that  these  global  variations 
preserve  the  shape  and  the  parameter  sets  (a, 0,0)  corresponding 
to  various  shapes  which  are  related  to  each  other  by  transforma¬ 
tions  like  scaling,  rotation,  translation  and  the  starting  point 


are  simply  related  by  Theorem  1.  On  the  other  hand,  the  para¬ 
meter  sets  of  quite  different  shapes  cluster  in  (a ,0,0)  far 
apart.  Fig.  3  illustrates  some  possible  transformations  of 
shapes.  The  constraint  on  the  parameter  set  (a, 6,0)  represent¬ 
ing  the  class  of  shapes  that  result  due  to  these  transformations 
is  given  in  Theorem  1. 

Theorem  1 :  Consider  a  closed  boundary  represented  by  r(l), 

r (2) , . . . , r (N)  ,  obeying  the  following  equation  involving  the 

coefficients  (a, 0,0): 

m 

r  (t )  =  a+  1  0  .  r  (t-t .  )  =/0u)(t),  t  =  1,2,.  .  .  ,N  (3.6) 

i=l  1  1 

All  other  boundaries  which  have  the  same  shape  as  the  boundary 
represented  by  r(l),  r (2) , . . . ,r (N)  can  be  generated  by  equation 
(3.6)  with  the  same  w ( • )  and  indexed  by  the  coefficients 
(a',0',0')  related  to  (a, 0,0)  by  (3.7): 

0'  =  0,  a'/ZF  =  a//0  (3.7) 

The  proof  is  given  in  the  appendix. 

The  significance  of  Theorem  1  is  as  follows:  Given  a  shape 
and  its  scaled  or  rotated  version,  the  corresponding  values  of 
0  and  a//0  are  identical.  Or  equivalently,  similar  shapes 
cluster  closely  in  the  (m+1)  parameter  space  (0 ,a//0  ).  Allowing 
variations  within  a  class  of  shapes  due  to  noise,  etc.,  we  see 
that  the  parameters  of  the  class  of  shapes  will  form  a  close 
cluster.  This  suggests  the  use  of  <J>  =  (0T,a//fT  )T  as  the 


feature  vector  for  classification  of  shapes.  This  will  be 
considered  in  detail  in  a  subsequent  paper. 


4.  Parameter  Estimation  and  the  Properties  of  the  Estimates 
The  unknown  parameters  (a, 6,6)  in  the  model  in  (3.2)  are 
estimated  from  the  observed  time  series  taken  from  a  given  boun¬ 
dary.  The  popular  method  of  estimation  of  parameters  is  the 
maximum  likelihood  estimate  (m.l.e).  The  m.l.e.  denoted  by 
(a*,0*,3*)  is  the  value  of  (a, 0,3)  which  maximizes  p(r(l),..., 
r(N)|a,0,3).  The  estimates  are  given  in  Theorem  2. 


Theorem  2 :  Let  (a*,0*,3*)  be  the  maximum  likelihood  estimates 
of  the  parameters  (a, 0,3).  Then 


N-i 

0*  =  Arg  | min  {-  £  in(l  -  0  *  ) 

k=0  ~  ~K 

0cR 


N  , 

+  N/2  in  E  (r  (t)  -  0  z(t  -  1)  -  ot(0))^} 
t=l 


(4.1) 


where 


«(0)  =  1/N  z  (r (t)  -  0Tz(t 
t=l  ~  ~ 


1)  ) 


a*  =  a (6*)  =  1/N  Z  (r(t)  -0*Tz (t  -  1)) 
t=l  ~  ~ 


(4.2) 


and 


N 

3*  =  1/N  Z  (r  (t) 
t=l 


a*  -  0*Tz(t  -  l))2 


(4.3) 


The  proof  is  given  in  the  appendix. 


8*  can  be  computed  by  a  gradient  algorithm.  For  reduced 
computations  we  suggest  the  use  of  the  least  square  (l.s.) 
estimates  given  below: 


(u'r,a)T 


N 

=  [  j:  u(t  - 

t=i~ 


l)uT(t 


-  1)1 


-1 


N 

T.  u  ( t  -  l)r(t) 
t=l~ 


(4.4) 


and 


K  =  (l/N)  >:  (r  (t)  -  (0T,a)u(t  -  l))2 
t=l 


(4.5) 


Also,  a  can  be  written  as 


a  *  l/N  l  (r (t)  -  6  z(t  -  1)) 
t=l  ~  ~ 


(4.6) 


where 

z (t  -  1)  =  (r(t  -  t.),...,r(t  -  t  ))T 
~  l  m 

The  l.s.  estimates  can  be  interpreted  as  obtained  by 
maximizing  p (r (1) , . . . ,r (N) | a ,9 , 6) ,  with  the  determinant  term 
| B ( 0 ) |  in  (3.4)  suppressed.  Since  the  relevant  matrix  B(0)  is 
nearly  lower  triangular  for  large  N,  the  l.s.  estimates  tend  to 
m.l.  estimates  asymptotically.  However,  the  easily  computable 
l.s.  estimates  also  possess  invariance  properties  w.r.t.  trans  - 
formations  like  rotation,  scaling  and  starting  point.  Numerical 


simulations  indicate  that  the  error  in  the  estimates  due  to 
the  suppression  of  the  determinant  term  is  negligible.  The 
details  of  the  simulations  are  given  in  Section  6. 

We  make  a  brief  note  regarding  the  estimation  of  parameters 
when  the  two-dimensional  representation  is  used  for  objects.  By 
making  the  simplifying  assumption  of  independence  of  the  set 
of  observations  (x^(*)}  and  (Xjt*))#  the  above  mentioned  methods 
of  estimation  are  directly  applicable.  More  details  on  the 
estimation  of  parameters  for  vector  processes  can  be  found  in  [12]. 

In  Theorem  1  we  considered  the  constraints  imposed  on  the 
parameter  sets  representing  the  class  of  shapes  that  result  due 
to  the  transformations  of  shape.  But  in  practice  we  work  with 
only  the  estimates  of  the  parameters.  It  is  then  necessary  to 
consider  the  constraints  imposed  on  the  estimates  of  the  para¬ 
meter  sets  from  the  shapes  that  result  due  to  transformations 
like  scaling,  rotation,  etc.  The  properties  possessed  by  the 
estimates  can  be  broadly  classified  as  those  related  to  the  shape 
and  those  related  to  the  statistics  of  the  data  derived  from  the 
shape.  For  instance,  the  invariance  of  (O,a//0)  under  transfor¬ 
mations  like  scaling  and  rotation  are  shape  related  properties, 
and  the  sufficiency  property  of  the  estimates  comes  under  the 
other  category.  The  estimate  (¥,'ci//B)  possesses  the  following 
properties : 

Pi)  Invariance  to  scaling  of  the  boundary. 

P2)  Invariance  to  variation  in  the  starting  point  used  in 

tracing  the  boundary. 


P3)  Invariance  to  rotation  of  the  boundary. 

P4)  Invariance  to  translation  of  the  boundary. 

P5)  The  l.s.  estimates  are  a  sufficient  statistic  for  the 
unknown  parameters. 

The  method  of  tracing  the  boundary  in  the  one-dimensional 
representation  is  such  that  an  effective  normalization  is  done 
w.r.t.  the  center  of  gravity  of  the  shape.  Thus  the  one-dimen¬ 
sional  representation  is  invariant  to  translation.  Since  the 
origin  of  the  two-dimensional  co-ordinate  system  is  fixed  at 
the  centroid  of  the  object  in  Method  1,  this  representation  is 
also  independent  of  translation. 

Also,  since  we  go  around  the  boundary  once  completely,  and 
use  the  entire  observation  set  in  estimating  the  parameters,  the 
effect  of  rotation  of  the  axes  and  the  effect  due  to  the  varia¬ 
tion  in  the  starting  point  are  the  same. 

Theorem  3:  Let  (a,  0,(5)  denote  the  estimates  of  the  Darameters  of 
the  model  representing  a  Type  1  object  and  let  (a^ and 
(a2'e2'®2*  represent  the  corresponding  estimates  of  the  model 
representing  a  non-Type  1  object.  Then  the  estimates  (9\a//6) 
and  {  ,  (02 ,a2//62)  }  satisfy  the  properties  P1-P4. 

The  proof  of  Theorem  3  is  given  in  the  appendix. 

We  state  a  theorem  below  regarding  P5.  A  proof  similar  to 
the  one  found  in  [13]  can  be  given  for  this  theorem. 


i 


Theorem  4:  The  l.s.  estimates  (a,9,B)  are  sufficient  statistics 
for  the  unknown  parameters  (a,9,B). 

There  are  two  consequences  resulting  from  the  existence  of 
sufficient  statistics  for  the  parameters  (a,9,B).  The  first 
consequence  is  that  all  the  information  contained  in  the  unknown 
parameters  (a,9,B)  is  also  contained  in  the  estimates  (a  .¥,3")  . 

The  second  consequence  is  that  (a, 9,6)  is  the  optimal 
feature  set  for  classification  purposes  [13].  Hence,  we  can 
restrict  our  attention  to  a  class  of  decision  rules  that  are 
functions  of  sufficient  statistics.  This  will  be  further  con¬ 
sidered  in  a  subsequent  paper. 

Choice  of  the  number  of  lag  terms;  We  have  so  far  done  the 
analysis  as  though  the  number  of  lag  terms  in  eq .  (3.2)  is  known. 
As  mentioned  in  Section  3,  we  need  not  have  all  the  lag  terms 
up  to  the  highest  lag  term,  the  lag  terms  present  being  indicated 
by  the  set  (t^L  The  problem  of  determining  the  number  of  lag 
terms  comes  under  a  general  class  of  problems  known  as  "system 
identification"  and  has  been  considered  elsewhere  [14].  Without 
going  into  the  details  we  simply  formulate  the  problem  and  give 
a  decision  rule  to  choose  the  best  m  for  the  given  boundary. 

Suppose  that  we  consider  three  possible  values  of  m,  ,  m2 , 
and  m3-  Corresponding  to  each  n^,  i  =  1,  2,  3,  we  will  write  the 
corresponding  circular  autoregressive  equation  E, ,  i  =  1,  2,  3: 


where  3,0  are  defined  in  (3.10)  and  (3.12)  and  (N)  is  some 
suitable  function  of  N.  Typically,  f i (N)  could  be  fcnN  (14),  or 
0.01.  The  decision  rule  to  choose  the  appropriate  m^  is  as  follows: 
Choose  the  model  k*  if 

k*  =  arg  min  (C,  }  (4.9) 

k  K 

The  choice  of  f ^ ( N )  in  (4.8)  is  a  tradeoff  between  the  two 
types  of  errors  that  arise  in  a  hypothesis  testing  problem. 

Consider  the  case  of  choosing  between  a  mfc^  and  an  (m+l)st  order 
model.  Then  it  can  be  shown  that  the  probability  of  choosing 
the  lower  order  model  when  indeed  the  higher  order  model  is  true 


(the  so-called  Type  II  error)  is  proportional  to  /f7(N).  When 
the  primary  interest  is  to  reconstruct  the  shape,  a  higher  order 
model  is  preferred  and  it  is  important  to  keep  the  Type  II  error 
under  control.  For  this  reason,  f^N)  should  be  0.01  or  smaller. 


f-. 

('■ 


5.  Coding  and  Reconstruction  Schemes 

In  this  section  we  discuss  the  coding  and  reconstruction 
schemes  possible  with  the  circular  model  developed  in  this 
paper.  We  consider  the  applicability  of  such  schemes  to  a 
class  of  two  tone  images  [2]  like  weather  maps  and  geographical 
contours,  typical  illustrations  of  which  are  shown  in  Fig.  4. 

The  given  closed  contour  is  represented  by  a  model  similar 
to  (3.2).  The  estimates  of  the  parameters  characterizing  the 
model  and  the  residuals  defined  as 


ui(t)  =  (1//  0))[r(t)  -  a  -  E  F.r(t  -  t  )] 

j=l  ]  3 


(5.1) 


are  computed.  The  CAR  model  considered  in  this  paper  splits 
the  information  contained  in  the  original  correlated  observations 
into  two  sets:  the  statistical  information  contained  in  the  set 
(a, 0,0)  and  the  information  contained  in  the  uncorrelated  resi¬ 
duals.  The  former  set  characterizes  a  class  of  shapes  with  dif¬ 
ferent  noise  sequences  while  the  combined  set  (ot ,  0",0  ,w  (1)  ,  . .  .  ,w(N) ) 
represents  a  particular  realization  of  the  shape.  In  actual 
implementation  the  exact  residuals  are  never  stored.  Either  they 
are  truncated  or  truncated  and  skipped,  i.e.,  some  of  the  trun¬ 
cated  residuals  are  replaced  by  those  from  a  (0,1)  Gaussian  random 
generator.  If  the  skipped  and  truncated  residuals  are  denoted 
by  { to  ( * )  >  and  (uiM*)}  respectively  and  { v  ( * )  >  denote  the  (0,1) 
Gaussian  random  numbers  then  one  method  of  generating  {£)(•)} 
would  be  as  in  (5.2): 


w ( 2 1 )  =  w*  (2t) , 


t  =  l,2,...,N/2 


u)(2t-l)  =  v  (2 1-1)  ,  t  -  l,3,...,N/2-l 

u)(2t+l)  =  w'  (2t+l)  ,  t  =  0,2,  .  . .  ,N/2-2 


(5.2) 


In  (5.2),  every  l/4th  of  truncated  residual  is  replaced  by  a 
(0,1)  Gaussian  random  number.  It  is  also  possible  to  generate 
{w(*)}  where  every  other  {<*>'(*)}  is  replaced  by  {v(*)K  The 
particular  scheme  is  generally  decided  by  the  degree  of  appro¬ 
ximation  desired.  Such  trimmed  residuals  are  quantized  and 
stored.  From  this  stored  information  a  closed  contour  that  is 
a  good  approximation  to  the  original  contour  can  be  generated 
by  the  simultaneous  scheme  described  below: 


Simultaneous  Scheme:  Given  N  trimmed  residuals  {£)(•)},  truncated 
or  truncated  and  skipped,  and  the  estimates  (a,F,jf),  find  r(i)  , 


i  *  1,2,. . . ,N  such  that  the  norm 
I  |B(0)R  -  A  -  $1  I  2 


(5.3) 


.  .  .  T  — T  —  —  — 

is  minimized,  where  R  =  (r (1) , . . . ,r (N) ) ,  A  =  [a, a,..., a],  Nxl 

A  /\  /V  T 

vector,  and  w  =  [u>  (1)  , . . .  ,o>  (N)  ]  .  An  optimal  solution  that  mini¬ 
mizes  the  expression  in  (5.3)  is  given  in  Theorem  5: 


Theorem  5 :  Let  u(l), (N)  be  the  trimmed  residuals  and 

(a, 0,6)  be  the  estimates  of  the  parameters  of  the  model.  Let 

k  =  Arg  {min( |  j  B ( 6 ) R  -  A  -  /6$| |^}}.  Then  $  can  be  computed 
R 

as  shown  below: 


A.,  % 

r  (k) 


N 


E  e^2lI/N>  ^  »-l)u,|i)(  k=lf2> 
i=l 


•  ,N,  j  -  /-I  (5.4) 


where 


u'  (i)  =  u(i)/X (i) 


(5.5) 


u  (i)  =  1/N  X  e‘j(27T/N)  (i  15  (n_1)  (a  +  /gfc(n>),  i=l,2,...,N  (5. 


n=l 


6) 


A(i)  =  (1  -  0>.)  ,  i  =  1,2,.  .  .  ,N 


(5.7) 


and 


4\.  =  (e 


-j  (2ir/N)  (i-Ut^  -j  (2u/N)  (i-1)  t2  -j  (2n/N)  (i-1)  t 


.  ,e 


m. 


The  proof  of  Theorem  5  is  given  in  the  appendix. 


Comments : 


1)  Reduction  in  the  computations  using  FFT  algorithms  is  possible 
due  to  the  following: 

a)  the  eigenvalues  of  the  circulant  matrices  can  be  written 
down  by  inspection  from  the  first  row  of  the  matrix,  and 

b)  The  eigenvectors  of  the  circulant  matrix  are  universal, 
i.e.,  all  possible  circulant  matrices  of  a  specified  order 
have  an  identical  set  of  Fourier  eigenvectors. 


a 


2)  Kq.  (5.7)  is  of  computational  complexity  0(N  log  N)  and 
so  is  Eq.  (5.9).  Since  other  computations  like  Eq.  (5.8)  are 
of  complexity  O(N),  the  overall  computational  complexity  is 
C)(N  log  N)  . 

Recursive  Scheme :  A  recursive  scheme  similar  to  the  DPCM 
scheme  can  be  suggested  for  the  reconstruction  of  contours. 

The  estimates  (a, 8, 8)  and  the  truncated  and  quantized  residuals 
are  used  in  (3.2)  along  with  m  initial  conditions  to  generate 
the  synthetic  observation  set.  Owing  to  the  uncorrelated  nature 
of  the  residuals  considerable  saving  in  bits  can  be  achieved. 

If  exact  residuals  are  used  the  original  closed  contour  can  be 
generated  but  in  practice  the  exact  residuals  are  never  stored. 
Only  truncated  and  quantized  residuals  are  stored.  Depending  on 
the  particular  truncation  and  quantization  schemes,  contours 
that  are  close  enough  to  the  original  contours  can  be  generated. 


6.  Results  of  Simulation 


In  this  section  we  summarize  the  results  of  simulations 
done  with  three  closed  boundaries  shown  in  Fig.  5,  denoted  as 
shapes  A,  B,  and  C.  Shapes  A  and  C  are  typical  Type  1  objects 
and  shape  B  is  a  typical  non-Type  1  object. 

6.1  Model  Order  Selection:  The  decision  rule  in  (4.9)  was  used 
to  choose  an  appropriate  model  for  different  values  of  f^(N). 

The  model  orders  chosen  for  different  values  of  f^(N)  for  shapes 
A  and  C  are  given  in  Table  1. 

Table  1.  Model  orders  for  different  f^(N)  in  (4.8). 


f.(N> 

MODEL  ORDER 

CHOSEN 

SHAPE  A 

SHAPE  C 

o.oi 

13 

15 

0.1 

13 

15 

1.0 

12 

6 

1.2 

6 

6 

1.5 

2 

2 

2.0 

2 

2 

As  discussed  before  by  decreasing  f^(N)  we  reduce  the  pro¬ 
bability  of  choosing  a  lower  order  model  when  the  higher  order 


model  is  true.  It  is  this  error  that  is  crucial  for  applications 
like  reconstruction.  The  value  of  f^  (N)  is  chosen  depending  on 
the  quality  of  reconstruction  desired. 

6.2  Accuracy  of  the  least  square  estimates:  To  check  the  accu¬ 
racy  of  the  l.s.  estimates  as  compared  to  m.l.  estimates,  the 
relative  contributions  of  the  first  term  and  the  quadratic  term 
in  (4.1)  to  the  log  likelihood  function  2,np  (r(l),,..,r(N)|a,0,B) 
were  computed  for  shapes  A  and  C,  for  a  twelfth  order  model,  and 
are  given  in  Table  1.  It  is  seen  that  the  contribution  of  log 
|det  B (0 )  |  is  negligible  compared  to  the  other  term. 

Table  2.  Relative  contributions  of  the  log  linear  term 
and  the  quadratic  term  in  (4.1). 


CONTRIBUTION 

CONTRIBUTION 

SHAPE 

DUE  TO  LOG  LINEAR 

DUE  TO  QUADRATIC 

TERN  IN  (4.1) 

TERM  IN  (4.1) 

A 

.1590 

24.96 

C 

-.2326 

48.90 

6.3  Reconstruction  of  boundaries:  We  give  examples  of  the  re¬ 
construction  of  boundaries.  Fig.  6  gives  the  results  of  simu¬ 
lation  experiments  with  shape  A.  Since  shape  A  is  a  Type  1 
object  the  one-dimensional  representation  was  used.  Least  square 


estimates  (a, 8,0)  were  computed.  In  Fig.  6  all  the  closed 
contours  have  been  reconstructed  using  the  simultaneous  scheme 
in  Theorem  5.  Fig.  6a  is  the  shape  reconstructed  using  a  second 
order  model  and  the  exact  residuals.  Fig.  6b  is  the  shape 
using  the  residuals  truncated  to  the  first  decimal  place.  The 
quality  of  reconstruction  is  close  to  Fig.  6a.  Figs.  6c  and  6d 
are  similar  to  Figs.  6a  and  6b  except  that  a  twelfth  order  model 
has  been  used.  All  the  shapes  in  Fig.  6  look  very  similar  to 
the  original  shape  A.  Then  the  natural  question  is  the  relevance 
of  the  decision  rule  used  in  fitting  an  appropriate  model.  This 
is  explained  as  follows:  the  residuals  of  an  inappropriate  model 
will  not  be  as  uncorrelated  as  those  of  an  appropriate  model 
and  hence  the  efficiency  of  coding  is  higher  when  a  model  of  the 
correct  order  is  used.  Since  the  data  generated  by  the  simulta¬ 
neous  model  is  circular  we  always  reconstruct  closed  contours. 
Fig.  7  gives  the  same  results  as  in  Fig.  6  for  shape  B. 

Fig.  8  gives  the  reconstructed  shape  A  using  the  recursive 
scheme  mentioned  earlier.  Fig.  8a  and  8b  are  the  results  of 
this  scheme  with  m  =  2  and  using  the  exact  and  the  truncated 
residuals  respectively.  Though  this  scheme  does  not  truly 
generate  a  closed  contour  (because  of  its  recursive  nature) 
still  the  quality  of  reconstruction  is  good.  Figs.  8c  and  8d 
are  similar  to  8a  and  8bf  with  a  twelfth  order  model.  Fig.  9 
gives  the  results  for  shape  B.  Fig.  9a  and  9c  are  reconstructed 
using  the  exact  residuals  and  m  =  2  and  12.  Fig.  9b  is  the 


result  of  the  recursive  scheme  with  truncated  residuals  and 


m  =  2.  Note  the  presence  of  a  spike,  due  to  the  noncircularity 

of  the  data.  Similarly,  Fig.  9d  is  the  reconstruction  with  m  =  12  j 

and  the  resulting  contour  is  not  closed.  ! 

As  mentioned  before,  the  circular  model  splits  the  informa-  ! 

tion  in  two  parts,  the  estimates  of  parameters  and  the  residuals.  j 

The  information  contained  in  the  estimates  represents  a  class  of 
shapes  corresponding  to  different  i.i.d  noise  sequences.  To  see  V 

how  much  information  is  contained  in  the  estimates,  synthetic  ^ 

generation  of  shapes  is  attempted.  The  results  are  given  in 

♦ 

Fig.  10  for  shape  A.  The  residuals  used  were  drawn  from  a  (0,1)  j_ 

Gaussian  random  number  generator.  Fig.  10a  and  10b  are  recon-  ! 

structed  with  m  =  2  and  the  recursive  and  simultaneous  schemes 

respectively.  Note  that  the  recursive  scheme  does  not  generate 

closed  contours.  Figs.  10c  and  lOd  are  for  m  =  12.  As  expected,  | 

the  quality  of  reconstruction  with  pseudo-random  numbers  is  not 

as  good  as  with  exact  residuals  but  the  general  shape  is  still 

recovered.  Also  note  that  models  of  order  12  generate  contours 

that  are  much  closer  to  the  original  than  the  contours  generated 

by  a  model  of  order  m  =  2. 

The  coding  schemes  mentioned  above  can  be  used  for  the  ! 

storage  of  contours  as  in  Fig.  4.  We. fit  a  stochastic  model 
to  every  closed  contour,  and  store  the  estimates  of  the  coeffi¬ 
cients  and  the  quantized  residuals.  For  contours  in  the  map 


that  are  not  closed,  we  can  fit  ordinary  autoregressive  models. 
The  map  can  be  reconstructed  by  reconstructing  each  contour 
separately. 


7.  Discussion 

In  this  section  we  compare  the  theory  developed  here  with 
other  well-known  methods  like  Fourier  analysis  [7-10,15].  The 
main  emphasis  of  the  theory  developed  here  is  to  demonstrate 
that  statistical  methods  can  be  developed  to  handle  a  wide  class 
of  closed  boundaries . 

In  the  familiar  Fourier  analysis  method  a  continuous  function 
is  extracted  from  the  boundary  and  represented  as  a  Fourier  series. 
These  methods  are  basically  deterministic  compared  to  the 
stochastic  models  considered  here.  Generally  stochastic  models 
are  found  to  be  more  useful  in  noisy  situations  and  tend  to 
describe  the  data  by  fewer  parameters  than  the  deterministic 
models  [16].  In  fact  by  using  the  theory  of  comparison  of  dynamic 
models  [14],  an  exact  comparison  of  the  Fourier  series  represen¬ 
tation  and  the  stochastic  model  considered  here  can  be  made. 

In  any  practical  situation,  the  Fourier  series  expansion  has 
to  be  truncated.  Because  of  the  truncation  in  the  Fourier  series 
representation,  the  data  reduction  from  the  observations  to  the 
Fourier  coefficients  is  not  exact.  We  have  given  a  decision  rule 
to  choose  the  appropriate  lag  terms  for  the  boundary  analysis. 

One  of  the  criticisms  of  current  boundary  coding  schemes 
[17,  p.  161]  is  that  the  points  (A,B)  and  (C,D)  in  Fig.  11 
which  are  geometrically  close  together  can  be  encoded  quite  far 
apart  in  the  resulting  string.  Such  boundaries  can  be  adequately 


represented  by  models  which  have  non-contiguous  lag  terms  [12]. 
Typically,  such  a  model  would  include  lag  terms  like  (x^ (t-1) , . . . , 
x^tt-N^)),  where  is  the  number  of  spaces  between  A  and  B. 

The  specific  CAR  model  considered  here,  because  of  its  cir¬ 
cular  nature,  is  intuitively  pleasing  for  representing  a  closed 
boundary.  The  two-dimensional  representation  considered  here  is 
used  to  represent  any  closed  curve  on  a  plane,  and  the  one-dimen¬ 
sional  representation  is  used  for  the  Type  1  objects.  The  time 
series  is  constructed  such  that  the  corresponding  shape  descrip¬ 
tion  is  invariant  to  the  transformations  of  the  boundary.  The 
vector  (0,a  //$  )  is  insensitive  to  these  transformations  of  the 
shape  and  so  is  a  natural  choice  for  shape  classification  prob¬ 
lems.  We  have  suggested  coding  and  reconstruction  algorithms. 

The  specific  simultaneous  reconstruction  algorithm  always  results 
in  a  closed  contour  in  O(NlogN)  operations.  We  have  also  con¬ 
sidered  the  problem  of  boundary  coding  for  archival  storage  of 
weather  maps  and  geographical  maps.  It  is  possible  to  extend 
these  methods  to  two-dimensional  images  [18]  where  the  edges  are 
coded  as  contours . 

Shape  descriptors  like  ratio  of  square  of  perimeter  to 
area  [19]  are  not  useful  for  synthesis  of  complicated  shapes. 

Also  they  are  not  robust  since  they  often  yield  similar  numerical 
values  for  contours  that  are  significantly  different  from  one 
another  [20]  ,  as  shown  in  Fig.  9. 


The  theory  developed  here  is  based  on  a  global  description 
of  shapes.  In  a  subsequent  paper  (Part  II)/  we  will  consider 
the  problems  of  shape  discrimination.  In  particular,  we  will 


be  interested  in  obtaining  decision  rules  that  are  invariant 
to  transformations  of  shapes. 


We  give  the  proofs  of  the  various  theorems  given  in  the 


body  of  the  paper. 

Theorem  1:  Given  a  closed  boundary  represented  by  (r(l),  r  (2 ),..., 
r(N)},  obeying  the  following  equation  involving  the  coefficients 
(a,0,B) , 

r  (t)  =  a  +  I  0 .  r  (t  -  t.)  +/Bu>(t),  t  =  1,2,... ,N  (1) 

i=l  1 

all  other  boundaries  which  have  the  same  shape  as  the  boundary 
represented  by  {r  (1) , . . . ,r (N) } ,  can  be  generated  by  equation  (1) 
with  the  same  w  (1)  , . . .  ,u  (N)  and  the  coefficients  (a',©'^')  re¬ 
lated  to  (a,6,B)  by  (2): 

0'  =  0,  =  a//ff  (2) 

Proof :  Shape  is  preserved  by  the  operations  of  scaling  and  rota¬ 
tion.  We  first  consider  the  constraints  imposed  on  the  parameters 
to  represent  a  class  of  shapes  that  differ  due  to  scaling, 
a)  Scaling ; 

Any  boundary  which  has  the  same  shape  as  {r(*)>,  but  is 
scaled,  is  generated  by  {kr(l),  kr  (2 ) , . . . ,kr  (N) } ,  where  k  is  the 
scale  factor,  and  obeys 


i 

.  j 


(3) 


m 

kr(t)  *  <*'  +  E  0!kr(t  -  t.)  +/P’,a)(t) 
i=l  1  1 


r  (t)  =  «'A  +  £  0!r(t  -  t .  )  +  (/fPyiOwft) 

i=l  1  1 


Since  r(t)  obeys  both  (1)  and  (4),  by  the  uniqueness  of 
solution  of  the  difference  equation,  (2)  follows, 
b)  Rotation:  Consider  the  rotated  boundary  r'(l 
such  that 


r' (i)  *  r (i+S) 


and 

<*>'  (i)  =  u>  (i+S) 


where  S  is  a  positive  integer.  Now 

m 

r' (t)  =  a*  +  I  0 !  r '  (t  -  t .  )  +  /3'u>’(t),  t  =  1,2,.. 
i=l  1  1 


or 


m 

r  (t  +  S)  =  a'  +  T.  0 !  r  (t  -  t.  +  S)  +  /(Putt  +  S)  ,  t  =  1 
i=l  1  1 


(4) 

the 

(N)  , 

(5) 

•  ,N  (6) 

,  2  ,  .  .  .  ,  N 


Equivalently , 


m 

r  (k)  =  a'  +  E  9 !  r  (k  -  t .  )  +  /fT'ojtk),  k  =  1+S,...,S+N  (7) 

i=l  1 

(1)  and  (7)  are  identical  and  hence  (2)  follows. 

Q.E.D. 

Proof  of  Theorem  2 :  From  (3.4),  taking  logarithms  and  using  (3.5), 
fcnp (r(l),...,r(N)|a,0,0)  =  -N/2|n2B 


-  1/26  E  (r  (t)  -  zT(t  -  1)9  -a)2  +  E  fln(l  -  ^0)  (8) 

t=l  ~  ~  k=0 


Differentiating  w.r.t.  a, 6  and  equating  to  zero, 


<x(0)  =  1/N  E  (r  (t)  -  zT(t  -  1)0) 
t=l 


0(0)  =  (1/N)  E  (r (t)  -  zT(t  -  1)0  -a ( 0 ) ) 2  (10) 

t=l  ~  ~ 


Substituting  (9)  and  (10)  in  (8)  we  have 


Unp (r(l),...,r(N) | a , 0 , 6)  »  -(N/2) | n (2n/N) 

N  7 

(-N/2)  £n  E  (r  (t )  -  zT(t  -  1)0  -  a(9)] 
t=i  ~  ~ 


+  E  fcn (1  -  Y*0)  -  N/2 
k=0  ~K~ 


Maximizing  in  p (r (1) , . . . , r (N) | a , 0 , 8)  given  above  w.r.t.  8  is 
equivalent  to  minimizing  (4.1)  w.r.t.  6.  Substituting  the  re¬ 
sulting  8*  in  (9)  and  (10)  we  obtain  Theorem  2. 

Theorem  3 :  Let  (a, 8, 8)  denote  the  estimates  of  the  parameters  of 
the  model  representing  a  Type  1  object  and  let  (a^,F.,F^)  and 
^*2'02'^2^  represent  the  corresponding  estimates  of  the  model 
representing  a  non  Type  1  object.  Then  the  parameter  sets 
( 0 ,  ct  / t/pj  and  {  (0^ ,a. /^T^)  ,  (62,02/^2))  satisfy  the  properties 
P1-P4. 

Proof:  (a)  Invariance  to  PI; 

Let 

r’ (t)  =  kr  (t)  ,  t  =  1,  2 ,  . .  .  ,N  (11) 

be  the  scaled  boundary. 

Let  the  l.s.  estimates  of  the  parameters  associates  with 
(r'  (•)>,  be  (a’ ,0* ,8*)  . 

From  (4.4)  and  (11),  it  follows  that  O'  =  0.  Also  from  (4.5) 
and  (4.6) 

(1/N) I  (r ' (t)  -  0Tz*(t  -  1)) 

a'/ZB^  =  3 - - — 1 - yj~2  ~  a/Z^  using  (11) 

1/N  l  (r'  (t)  -  e^z'tt  -  1)  -  a(9))2 
t=l 

(b)  Invariance  to  P2 :  for  this  case  r' (t)  =  r(t+S)  where  S  is 
the  unit  of  rotation.  By  actually  substituting  in  the  expressions 
for  O'  and  and  using 


i 


r (N+k)  =  r(k), 

it  is  shown  that  0 '  =  F  and  a'/ZB'  =  a/ZF. 

Proof  of  Theorem  5; 

From  (5.3)  we  have 

&  =  B(0)-1[A  +/S  (1.2) 

where  8(F)  is  a  NxN  circulant  matrix.  It  is  well  known  [11] 
that  the  eigenvectors  of  a  NxN  circulant  matrix  are  the  Fourier 
vectors  >  where  4^ ,  the  jth  component  of  the  ith 

eigenvector  is  given  by 

_  /=T(2u/N) (i-1) (j-1)  (13) 

ij 

and  the  distinct  eigenvalues  A (1) , . . . , A (N)  are  given  by 


A  (k)  =  (1  -  0T^k) 


(1*) 


b  =  tA  +  /gw] 

Expanding  b  in  terms  of  the  Fourier  vectors  { 4>  we  have 

N 

b  =  Z  a.  <p  . 

'  i=l  1~1 


(15) 


(le: 


where 


«i  =  (J^b/N 


(27) 


‘J 

t"i 


and  denotes  the  complex  conjugate  of  <|>T.  Hence  from  (12), 
(15),  (16),  and  (17) 

N 

&  *  B(0)-1b  =  B(F)"1  Z  .  <j>*b/N  (18 


Using 

B(e)'1^  =  ( 1/ A  (i) )  4>  (i) 


the  LHS  of  (18)  is  equal  to 

N 

=  (19) 

Using  (13),  (14),  and  (19),  the  k*"*1  component  of  $  is 

A  N 

r  (k)  =  Z  exp  [  i/-12tt  (i  -  1)  (k  -  1)/N]/(1  -  S_T'l' .  ) 
i=l  ~  -1 


(  Z  exp [-  /-12ir  (i  -  1)  (n  -  1)/N](a  +/lw(n))  (20) 

n=l 

Defining 

N 

u(i)  »  1/N  Z  exp  (-  /=12ir(i  -  1)  (n  -  1)/N)  (a  +  /l^j(n)) 
n=l 


and 

u’  (i)  =  u  (i)  /  (1  -  0T'Fi ) 

we  have 

A  N 

t(k)  =  Z  exp[/^l2ir  (i-1)  (k-l)/N]u’  (i)  ,  k  =  1,2,.  . .  ,N  (23) 
i=l 


Q  .  E  .  D  . 
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Fig.  3.  Variations  In  shape:  a)  original  shape,  b)  noisy  shape  (local  trans¬ 
formations),  c)  original  shape  globally  transformed  (scaling). 


Fig.  4.  Contour  maps. 


Fig.  6.  Reconstruction  of  shape  A  using  the  simultaneous  scheme,  a)  m-2 
vlth  exact  residuals,  b)  m*2  with  residuals  truncated  to  first 
decimal  place,  c)  similar  to  a)  with  m-12,  d)  similar  to  b) 
with  m*12.  Note  that  b)  and  d)  are  close  to  the  original. 


Fig.  8.  Reconstruction  of  shape  A  using  the  recursive  scheme,  a)  m*2  with 
exact  residuals,  b)  m-2  with  residuals  truncated  to  first  decimal 
place,  c)  similar  to  a)  with  m-12,  d)  similar  to  b)  with  m-12. 


Fig.  10.  Synthetic  generation  of  shape  A  with  residuals  drawn  from  a  (0  1) 
pseudorandom  Gaussian  generator,  a)  m«2  with  recursive  scheme, 
b)  ro-2  with  simultaneous  scheme,  c)  similar  to  a)  with  m-12, 
d)  similar  to  b)  with  m“12.  Note  that  the  contours  constructed 
with  the  recursive  scheme  are  not  closed. 
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